# Wczytanie bibliotek
from pathlib import Path
from time import strftime
import pandas as pd
import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
import seaborn as sns
import plotly.express as px
import tensorflow as tf
import keras_tuner as kt
import country_converter as cc
from pandas.plotting import scatter_matrix
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler, OneHotEncoder
from sklearn.base import BaseEstimator, TransformerMixin
from sklearn.model_selection import train_test_split, GridSearchCV
from sklearn.ensemble import RandomForestRegressor
from sklearn.tree import DecisionTreeRegressor
from sklearn.linear_model import LinearRegression
from sklearn.svm import SVR
from sklearn.metrics import mean_squared_error, r2_score, mean_absolute_error
from xgboost import XGBRegressor
import shap
import joblib
import warnings
from classes_preprocessing import (ColumnDropper, ChangeNameCountry, AddRegion,TransPower, TransLog, OneH, ChangeToBinnary)

1. Cel projektu

Głównym celem badania jest opracowanie modelu uczenia maszynowego służącego do predykcji poziomu szczęścia w poszczególnych krajach na podstawie wskaźników ekonomicznych, społecznych oraz regionalnych. Dodatkowy cel stanowi weryfikacja postawionych w toku pracy pytań oraz hipotez badawczych.

2. Charakterystyka zbioru danych

Źródło danych

Dane zostały zaczerpnięte z platformy kaggle. Wykorzystany zbiór zawiera zestawienie wyników z lat 2015–2024, dla większości państw świata. Wyjątek stanowią kraje, dla których dane nie były dostępne, są to:

• Andora, Antigua i Barbuda, Dominika, Kiribati, Korea Północna, Liechtenstein, Malediwy, Wyspy Marshalla, Monako, Nauru, Palau, Saint Kitts i Nevis, Saint Lucia, Saint Vincent i Grenadyny, San Marino, Sudan, Tonga, Tuvalu, Watykan

Oryginalnie dane pochodzą ze Światowego Raportu Szczęścia (World Happines Report) - oficjalnego, globalnego, rankingu szczęścia, publikowanego corocznie przy wsparciu ONZ. Raport ten jest uznawany za najważniejsze źródło wiedzy o jakości życia i dobrostanie społecznym na świecie.

Opis danych

Zmienna	Znaczenie
Ranking	Zmienna porządkowa wskazująca pozycję kraju w globalnym zestawieniu szczęścia
Happiness score	Wartość szczęścia [wyznaczana na podstawie ankiet]
GDP per capita	Produkt krajowy brutto na mieszkańca skorygowany o parytet siły nabywczej
Social support	Poziom wsparcia społecznego [średnia z odpowiedzi binarnych w ankietach]
Healthy life expectancy	Średnia oczekiwana długość życia w zdrowiu
Freedom to make life choices	Postrzegana swoboda podejmowania ważnych decyzji życiowych [procent osób zadowolonych wg ankiet]
Generosity	Poziom hojności w oparciu o ostatnie darowizny [na podstawie ankiet, skorygowane o PKB]
Perceptions of corruption	Postrzegana korupcja w sektorze rządowym i biznesowym [średnia z odpowiedzi binarnych w ankietach]
Year	Rok obserwacji
Country	Nazwa kraju
Regional indicator	Region geograficzny

3. Czyszczenie danych

# Wczytanie danych pochodzących z platformy Kaggle
data = pd.read_csv("world_happiness_combined.csv", sep=";")

Nadanie odpowiednich typów

# Wyświetlenie informacji o ilości i typie każdej zmiennej
data.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 1502 entries, 0 to 1501
Data columns (total 11 columns):
 #   Column                        Non-Null Count  Dtype 
---  ------                        --------------  ----- 
 0   Ranking                       1502 non-null   int64 
 1   Country                       1502 non-null   object
 2   Regional indicator            1499 non-null   object
 3   Happiness score               1502 non-null   object
 4   GDP per capita                1502 non-null   object
 5   Social support                1502 non-null   object
 6   Healthy life expectancy       1502 non-null   int64 
 7   Freedom to make life choices  1502 non-null   object
 8   Generosity                    1502 non-null   object
 9   Perceptions of corruption     1502 non-null   object
 10  Year                          1502 non-null   int64 
dtypes: int64(3), object(8)
memory usage: 129.2+ KB

Zmienne ilościowe zostały wczytane jako typ object z powodu występowania przecinków w danych. Z tego względu wykonano konwersje niektórych zmiennych na typ float.

# Lista kolumn którym należy zmienić typ
fix_cols = [
    'Happiness score', 
    'GDP per capita', 
    'Social support', 
    'Freedom to make life choices', 
    'Generosity', 
    'Perceptions of corruption'
]
# Zmina typów kolumnom z listy
for col in fix_cols:
    data[col] = data[col].astype(str).str.replace(',', '.')
    data[col] = pd.to_numeric(data[col], errors='coerce')

# Wyświetlenie typów po zmianach
data.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 1502 entries, 0 to 1501
Data columns (total 11 columns):
 #   Column                        Non-Null Count  Dtype  
---  ------                        --------------  -----  
 0   Ranking                       1502 non-null   int64  
 1   Country                       1502 non-null   object 
 2   Regional indicator            1499 non-null   object 
 3   Happiness score               1502 non-null   float64
 4   GDP per capita                1502 non-null   float64
 5   Social support                1502 non-null   float64
 6   Healthy life expectancy       1502 non-null   int64  
 7   Freedom to make life choices  1502 non-null   float64
 8   Generosity                    1502 non-null   float64
 9   Perceptions of corruption     1502 non-null   float64
 10  Year                          1502 non-null   int64  
dtypes: float64(6), int64(3), object(2)
memory usage: 129.2+ KB

Oczyszczenie danych z błędów grubych

Podczas eksploracyjnej analizy danych zidentyfikowano błędy grube występujące w latach 2015-2018 w kolumnie GDP per Capita.

Poniższy wykres na przykładzie Bangladeszu ilustruję skale problemu. W rzeczywistości w 2016 roku nie odnotowano w tym kraju żadnych gwałtownych zmian ekonomicznych, które uzasadniałyby tak nagły skok, a nastepnie spadek wartości. Podczas analizy natknięto się na podobne błędy w innych krajach w latach 2015-2018. Stąd też podjęto decyzję o ograniczeniu analizowanego zbioru do lat 2019-2024, gdyż są to dane najlepszej jakości.

# Wybranie danych dla Bangladeshu
data_Bnagladesh = data[data['Country'] == 'Bangladesh'].sort_values('Year')
# Tworzenie wykresu
sns.lineplot(data=data_Bnagladesh, x='Year', y='GDP per capita', color='gray', marker='o')
point_2016 = data_Bnagladesh[data_Bnagladesh['Year'] == 2016]
sns.scatterplot(data=point_2016, x='Year', y='GDP per capita', color='red', s=100, zorder=5)
plt.title('Zmiana PKB per capita w Bangladeszu')
plt.show()

# Odrzucenie obserwacji przed 2019 rokiem
data = data[data['Year'] >= 2019]

Zidentyfikowano także błąd w nazwie kraju: zmieniono 'Argelia' (zapis hiszpański) na 'Algeria' (zapis angielski), dostosowując rekord do reszty zbioru.

# Zmiana nazwy kraju 'Argelia' na 'Algeria'
data['Country'] = data['Country'].replace('Argelia', 'Algeria')

Ujednolicenia klasyfikacji regionalnej

Podczas analizy zmiennych kategorycznych wykryto niespójności w kolumnie Region indicator. Zauważono przypadki, w których ten sam kraj w zależności od roku pomiaru zostawał przypisany do różnych regionów. W tej kolumnie występowały także 3 braki.

# Wyświetlenia danych dla Cypru
print(data.loc[data['Country'] == 'Cyprus', ['Country', 'Year', 'Regional indicator']])

     Country  Year  Regional indicator
673   Cyprus  2019                 NaN
824   Cyprus  2020      Western Europe
963   Cyprus  2021      Western Europe
1110  Cyprus  2022  Sub-Saharan Africa
1270  Cyprus  2023  Sub-Saharan Africa
1391  Cyprus  2024      Western Europe

W celu naprawy błędów przypisań regionów, wykorzystano bibliotekę country_conveter do wygenerowanie poprawnych i jednolitych dla każdego kraju regionów.

# Dodanie krótkiej nazwy kraju (To na niej będzie opierało się przypisanie regionu)
data['Country_cc'] = cc.convert(names=data['Country'], to='name_short', not_found='not_found')
# Dodanie regionu
data['Region_cc'] = cc.convert(names=data['Country_cc'], to='UNregion')
# Odrzucenie państw które bazowo nie występowały w zbiorze
data = data[~data['Country'].isin(['North Cyprus', 'Somaliland region'])]
# Usunięcie starych kolumn
data = data.drop(columns=['Country', 'Regional indicator'])
# Usunięcie "_cc" z nazw nowych kolumn
data.columns = data.columns.str.replace('_cc', '', regex=False)

Odrzucenie zbędnych kolumn

W ramach czyszczenia danych odrzucono kolumnę 'Ranking'. Zmienna ta powiela informacje z kolumny 'Happines score' (bezpośrednio z niej powstaje), oraz nie stanowi cechy kraju.

# Usunięie kolumny Ranking
data = data.drop(columns=['Ranking'])

4. Pytania badawcze i hipotezy

Pytania:

• Czy możliwe jest na podstawie zaledwie kilku zmiennych określić poziom szczęścia w danym państwie?

• Czy pandemia COVID-19 negatywnie wpłynęła na poziom sczęścia?

Hipotezy:

• Zamożność społeczeństwa, wyrażona przez PKB na mieszkańca skorygowana o parytet siły nabywczej, stanowi jeden z najsilniejszy predyktorów poziomu szczęścia

• Hojność jest najsłabszym predyktorem szczęścia

5. Eksploracyjna analiza danych (EDA)

Analiza korelacji

# Wybranie tylko kolumn numerycznych
numeric = data.select_dtypes(include=['number'])
# Obliczenie korelacji
corr_matrix = numeric.corr()
# Tworzenie wykresu korelacji Pearsona
plt.figure(figsize=(8, 6))
sns.heatmap(corr_matrix, 
            annot=True,
            fmt=".2f",
            cmap='coolwarm',
            vmin=-1, vmax=1,
            linewidths=0.5,
            square=True)
plt.title("Macierz Korelacji Pearsona", fontsize=16)
plt.tight_layout()
plt.show()

Na podstawie powyższego wykresu można wyciągnąć następujące wnioski:

• Główne filary szczęścia: Analiza wykazała, że GDP per capita oraz Social support mają najsilniejszy, pozytywny związek z Happiness score. Sugeruje to, że zamożność i poczucie bezpieczeństwa w społeczeństwie to kluczowe elementy, które pozwalają najlepiej przewidzieć poziom szczęścia w danym kraju.

• Rola zdrowia i wolności: Istotny wpływ na szczęście mają również Healthy life expectancy oraz Freedom to make life choices. Są one wyraźnie skorelowane liniowo z Happiness score, choć nieco słabiej niż czynniki ekonomiczne.

• Brak istotnej zależności liniowej dla hojności i korupcji: Zmienne takie jak Generosity czy Perceptions of coruption nie wykazują wyraźnej liniowej zależności z Happiness score. Może to oznaczać, że ich wpływ jest bardziej złożony lub zależy od innych czynników, czego prosta analiza korelacji nie wychwyciła.

• Współliniowość: Zauważalna jest silna zależność między zmiennymi objaśniającymi – na przykład zmienna GDP per capita jest silnie liniowo skorelowana ze zmiennymi Social support, oraz Healthy life expectancy. Sugeruje to istotną zależność: kraje w których ludzie są bogatsi charakteryzują się zazwyczaj lepszym wsparciem społecznym oraz wyższą oczekiwaną długością życia obywateli. Taka wpółliniowość może jednak w przyszłości potencjalnie źle oddziaływać na niektóre modele uczenia maszynowego.

• Ukryty koszt korupcji: Ciekawym odkryciem jest wyraźna ujemna korelacja liniowa między Healthy life expectation a Perceptions of corruptions. Mimo że korupcja nie wpływa bezpośrednio na samo poczucie szczęścia w tym zestawieniu, to dane sugerują, że w państwach o wyższym wskaźniku korupcji ludzie żyją krócej.

Analiza trendu

# Przygotowanie dodatkowej ramki pomagającej zaznaczyć lata pandemii
covid_data = data[data['Year'].isin([2020, 2021, 2022])]
covid_means = covid_data.groupby('Year')['Happiness score'].mean().reset_index()
# Tworzenie wykresu
plt.figure(figsize=(10, 6))
sns.lineplot(data=data, x='Year', y='Happiness score', 
             marker='o', errorbar=None, color='gray', label='Średni wynik')
sns.scatterplot(data=covid_means, x='Year', y='Happiness score', 
                color='red', s=150, zorder=2, label='Lata największych obostrzeń związanych z pandemią (2020-2022)')
plt.title('Globalny trend poziomu szczęścia (wyróżnienie lat pandemii)')
plt.ylabel('Średni wynik (Happiness Score)')
plt.xlabel('Rok')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

Analiza trendu globalnego poziomu szczęścia w latach 2019–2024 wskazuje na stabilność tego wskaźnika. Średnia wartość Happiness Score oscyluje w wąskim przedziale (ok. 5,40 – 5,56), co sugeruje, że w skali świata deklarowany poziom szczęścia nie uległ drastycznym zmianom.

Szczególną uwagę zwraca okres pandemii COVID-19. Wbrew intuicyjnym przypuszczeniom, w latach największego nasilenia restrykcji (2020–2022) nie odnotowano załamania trendu. Co więcej, widoczny jest wręcz minimalny wzrost średniego poczucia szczęścia względem roku bazowego 2019. Pozwala to na sformułowanie odpowiedzi na postawione wcześniej pytanie badawcze: pandemia COVID-19 nie miała istotnego, negatywnego wpływu na globalny poziom szczęścia.

Potwierdzenie braku istonych różnic między latami

Za pomocą testu ANOVA można potwierdzić wyciągniętą na podstawie wcześniejszgo wykresu tezę o braku istotnych różnic w Happiness score między latami.

Weryfikacja założeń:

Założenie o jednorodności wariancji: W celu weryfikacji, czy rozrzut wyników jest zbliżony we wszystkich badanych latach, przeprowadzono test Levene'a. Sformułowano następujące hipotezy:

• $H_0$: Wariancje zmiennej Happiness Score są równe we wszystkich latach
• $H_1$: Przynajmniej jedna z wariancji różni się istotnie od pozostałych

# Przygotowanie danych
years = sorted(data['Year'].unique())
groups = [data[data['Year'] == year]['Happiness score'] for year in years]
# Przeprowadzenie testu Levene'a
F, p = stats.levene(*groups)
# Wyświetlenie p_wartośći z testu
print(f"P-wartość z przeprowadzonego testu Levene'a wynosi: {p:.3f}")

P-wartość z przeprowadzonego testu Levene'a wynosi: 0.944

Na każdym sensownym przyjętym poziomie istotności nie ma podstaw do odrzucenia $H_0$.
Oznacza to, że założenie o jednorodności wariancji zostało spełnione.

Założenie o normalności rozkładów danych: W celu weryfikacji, czy zmienna Happiness Score posiada rozkład normalny w poszczególnych latach, przeprowadzono test Shapiro-Wilka. Badanie wykonano niezależnie dla każdego roku. Przyjęto następujące hipotezy:

• $H_0$: Rozkład zmiennej Happiness Score w analizowanym roku jest zgodny z rozkładem normalnym.
• $H_1$: Rozkład zmiennej Happiness Score w analizowanym roku różni się istotnie od rozkładu normalnego.

# Przeprowadzenie testu Shapiro-Wilka dla każdego roku
for year, group_data in zip(years, groups):
    S, p = stats.shapiro(group_data)
    print(f"Rok {year}: p-wartość = {p:.3f} ")

Rok 2019: p-wartość = 0.150 
Rok 2020: p-wartość = 0.256 
Rok 2021: p-wartość = 0.458 
Rok 2022: p-wartość = 0.365 
Rok 2023: p-wartość = 0.077 
Rok 2024: p-wartość = 0.004 

Z powyższych wyników testów Shapiro-Wilka można wysnuć następujące wnioski:
Przyjmując standardowy poziom istotności $\alpha$ = 0,05 wszystkie lata oprócz roku 2024 spełniają założenie o normalności rozkładu zmiennej Happiness score.
Mimo tego odstępstwa, biorąc pod uwagę dużą liczebność próby oraz wynik testu Levene'a, który wykazał jednorodność wariancji w grupach, zdecydowano się przeprowadzić parametryczny test ANOVA.

# Wyświetlenie ilości obserwacji w każdym roku
(data['Year'].value_counts().sort_index().to_frame(name='Liczba obserwacji').T.style.set_caption("Liczebność próby w poszczególnych latach").set_table_styles([{'selector': 'caption', 'props': [('font-weight', 'bold'), ('font-size', '14px')]}]))

Liczebność próby w poszczególnych latach

Year	2019	2020	2021	2022	2023	2024
Liczba obserwacji	155	152	148	145	137	140

Weryfikacja tezy o braku istotnych różnic zmiennej Happiness score między latami : W celu weyfikacji tezy przeprowadzono test ANOVA. Sformułowano nastęujące hipotezy:

• $H_0$: Średnie wartości zmiennej Happiness score są równe we wszystkich latach
• $H_1$: Średnia wartość Happiness Score w co najmniej jednym roku różni się istotnie od pozostałych.

# Przeprowadzenie testu ANOVA
F, p = stats.f_oneway(*groups)
# Wyświetlenie p-wartości
print(f"P-wartość testu ANOVA wynosi: {p:.3f}")

P-wartość testu ANOVA wynosi: 0.858

Na żadnym sensownym przyjętym poziomie istotności nie ma podstaw do odrzucenia $H_0$.
Można zatem przyjąć, że średnie we wszystkich latach są sobie równe. Test ten potwierdza, że pandemia COVID-19 nie miała istotnego wpływu na poziom szczęścia.
Z tego względu zdecydowano się na nieuwzględnianie zmiennej Year w modelach uczenia maszynowego. Wprowadzałaby ona jedynie niepotrzebny szum, a dzięki jej pominięciu stworzone modele staną się bardziej uniwersalne.

Analiza regionalna

# Wyświetlanie tabeli z ilośćią obserwacji dla każdego regionu
(data['Region'].value_counts().to_frame().T.style.set_caption("Liczba obserwacji w każdym regionie").set_table_styles([{'selector': 'caption', 'props': [('font-weight', 'bold')]}]))

Liczba obserwacji w każdym regionie

Region	Western Asia	Western Africa	Southern Europe	Eastern Africa	South America	Northern Europe	Eastern Europe	South-eastern Asia	Southern Asia	Central America	Western Europe	Eastern Asia	Middle Africa	Northern Africa	Southern Africa	Central Asia	Caribbean	Australia and New Zealand	Northern America
count	87	84	78	76	60	60	58	54	45	42	42	36	30	29	28	27	17	12	12

# Tworzenie słownika aby przypisać odpowiedni kontynent na podstawie regionu
region_to_continent = {
    'Northern Europe': 'Europe',
    'Western Europe': 'Europe',
    'Australia and New Zealand': 'Oceania',
    'Northern America': 'America',
    'Central America': 'America',
    'Western Asia': 'Asia',
    'Eastern Europe': 'Europe',
    'Southern Europe': 'Europe',
    'Eastern Asia': 'Asia',
    'South America': 'America',
    'South-eastern Asia': 'Asia',
    'Caribbean': 'America',
    'Central Asia': 'Asia',
    'Eastern Africa': 'Africa',
    'Southern Asia': 'Asia',
    'Northern Africa': 'Africa',
    'Western Africa': 'Africa',
    'Middle Africa': 'Africa',
    'Southern Africa': 'Africa'}
# Tworzenie nowej ramki ze zmienna kontynent
data_continent = data.copy()
data_continent['Continent'] = data_continent['Region'].map(region_to_continent)
# Tworzneie wykresy boxplot względem szcęścia z podziałem na kontynenty i regiony
plt.figure(figsize=(12, 6))
sns.boxplot(data=data_continent, 
            x='Region', 
            y='Happiness score',
            hue='Continent',
            dodge=False)
plt.xticks(rotation=45, ha='right')
plt.title('Rozkład poziomu szczęścia w poszczególnych regionach', fontsize=16)
plt.xlabel('Region')
plt.ylabel('Happiness Score')
plt.legend(title='Continent', bbox_to_anchor=(1.05, 1), loc='upper left')
plt.tight_layout()
plt.show()

Analiza rozkładu zmiennej Happiness score w podziale na regiony pozwala na wyodrębnienie 3 wyraźnich grup:

• Liderzy rankingu (Wysoki poziom szczęścia): Do tej grupy należą zaledwie cztery regiony: Europa Północna, Europa Zachodnia, Ameryka Północna oraz Australia i Nowa Zelandia. Cechują się one zdecydowanie najwyższą medianą poziomu szczęścia w porównaniu do reszty świata
• Grupa umiarkowana (Średni poziom szczęścia): Jest to najliczniejsza i najbardziej zróżnicowana grupa, obejmująca pozostałe regiony Europy (Środkowo-Wschodniej), Ameryki (Centralną oraz Południową) oraz większość Azji. Regiony te plasują się pośrodku, prezentując przeciętne wartości wskaźnika szczęścia.
• Regiony o najniższych wskaźnikach: Trzecią grupę stanowią obszary, w których odnotowano najniższy poziom szczęścia. Zaliczają się do niej wszystkie regiony Afryki oraz Azja Południowa.

# Tworzenie interaktywnej mapy
fig = px.choropleth(data, locations="Country", locationmode='country names', color="Happiness score", hover_name="Country", color_continuous_scale=px.colors.sequential.Plasma, title="Mapa poziomu szczęścia na świecie")
fig.update_layout(height=600,width=1000,)
fig.show()

Interaktywna mapa stanowi wizualne dopełnienie wcześniejszej analizy, potwierdzając wyraźne dysproporcje w poziomie szczęścia pomiędzy poszczególnymi regionami świata.

Anliza rozkładów zmiennych

# Wyświetlenia hisotgramów rozkładów zmiennych
data.hist()
plt.tight_layout(pad=1.0)
plt.show()

Analiza histogramów wykazała istotną skośność rozkładów wybranych zmiennych, co może negatywnie wpłynąć na wyniki modeli uczenia maszynowego. W celu zbliżenia rozkładów do normalnego zastosowano następujące transformacje (w ramach pipeline modeli):

• Dla rozkładów lewoskośnych (Social Support, Freedom to make life choices, GDP per Capita) zmienne podniesiono do kwadratu
• Dla rozkładów prawoskośnych (Generosity) zmienne spierwiastkowano

Analiza rozkładu zmiennej Perceptions of corruption pokazuje, że społeczeństwo jest zazwyczaj zgodne, co do występowania korupcji w swoim państwie – jest ona postrzegana albo jako zjawisko powszechne, albo marginalne.

# Wyliczanie jaka ilość obserwacji mięści sie w przedziale: <0, 0.3> u <0.7, 1>
ilosc = ((data['Perceptions of corruption'] <= 0.3) | (data['Perceptions of corruption'] >= 0.7)).sum()
# Wyliczanie jaki to % całego zbioru
procent = (ilosc / len(data)) * 100
# Wyświetlenie wyniku
print(f"Potwierdzają to dane: aż {procent:.2f}% obserwacji znajduje się na krańcach skali (w przedziale <0, 0.3> u <0.7, 1>).")

Potwierdzają to dane: aż 77.77% obserwacji znajduje się na krańcach skali (w przedziale <0, 0.3> u <0.7, 1>).

W związku z tym, aby nie wprowadzać zbędnego szumu do modeli uczenia maszynowego poprzez wartości pośrednie, zdecydowano się na przekształcenie tej cechy w zmienną binarną

Rozkłady po zmianie:

# Tworzenie nowej ramki aby pokazać jak zmieni się rozkład nie modyfikując oryginalnej ramki
data_transformed = data.copy()
# Zmiana rozkładów
data_transformed['Social support'] = data_transformed['Social support'] ** 2
data_transformed['Freedom to make life choices'] = data_transformed['Freedom to make life choices'] ** 2
data_transformed['Generosity'] = np.sqrt(data_transformed['Generosity'])
data_transformed['GDP per capita'] = data_transformed['GDP per capita'] ** 2
data_transformed['Perceptions of corruption'] = data_transformed['Perceptions of corruption'].apply(lambda x: 0 if x < 0.5 else 1)
# Wyświetlenie histogramów rozkładów zmiennych
data_transformed.hist()
plt.tight_layout(pad=1.0)
plt.show()

6. Tworzenie modeli uczenia maszynowego

# Wczytanie oryginalnego zbioru danych. 
# Wszelkie transformacje (jak np. zmiana rozkładów) zostaną wykonane wewnątrz Pipeline'u 
# Takie podejście gwarantuje automatyczne przetworzenie przyszłych, nowych danych.
data = pd.read_csv("world_happiness_combined.csv", sep=";")

# Wstępne czyszczenie danych
# Poniższe operacje naprawiają błędy specyficzne dla tego konkretnego pliku i nie są częścią głównego Pipeline'u predykcyjnego
# ponieważ zakładamy, że nowe dane będą poprawne.

# Ograniczenie sie do danych od 2019 roku, gdyż w oryginalnych danych wcześniejsze obserwacje zawierały błędy
data = data[data['Year'] >= 2019]

# Zmiana nazwy kraju 'Argelia' na 'Algeria'
data['Country'] = data['Country'].replace('Argelia', 'Algeria')

# Zmiana typów niektórych kolumn
# Lista kolumn którym należy zmienić typ
fix_cols = [
    'Happiness score', 
    'GDP per capita', 
    'Social support', 
    'Freedom to make life choices', 
    'Generosity', 
    'Perceptions of corruption'
]
# Zmina typów kolumnom z listy
for col in fix_cols:
    data[col] = data[col].astype(str).str.replace(',', '.')
    data[col] = pd.to_numeric(data[col], errors='coerce')

# Oddzielenie zmiennej objaśniającej od zmiennych objaśnianych
X = data.drop('Happiness score', axis=1)
y = data['Happiness score']
# Dzielenia zbioru na zbiór testowy i treningowy
X_train_1, X_test_1, y_train_1, y_test_1 = train_test_split(X, y, test_size=0.2)

# Inicjalizacja listy do przechowywania wyników wszystkich modeli
model_results = []
# Funkcja pomocnicza do obliczania i zapisywania metryk modelu
def add_model_summary(model_name, y_train_true, y_train_pred, y_test_true, y_test_pred):
    # Obliczanie metryk dla zbioru treningowego
    train_r2 = r2_score(y_train_true, y_train_pred)
    train_rmse = np.sqrt(mean_squared_error(y_train_true, y_train_pred))
    train_mae = mean_absolute_error(y_train_true, y_train_pred)
    # Obliczanie metryk dla zbioru testowego
    test_r2 = r2_score(y_test_true, y_test_pred)
    test_rmse = np.sqrt(mean_squared_error(y_test_true, y_test_pred))
    test_mae = mean_absolute_error(y_test_true, y_test_pred)
    # Zapisanie wyników do słownika
    result = {
        "Model": model_name,
        "Train R2": train_r2,
        "Test R2": test_r2,
        "Train RMSE": train_rmse,
        "Test RMSE": test_rmse,
        "Train MAE": train_mae,
        "Test MAE": test_mae
    }
    # Dodanie danych do listy
    model_results.append(result)

Model regresyji liniowej

# Tworzenie pipeline do regresji liniowej [Wykorzystano autorskie transformery zdefiniowane w zewnętrznym pliku "classes_preprocessing.py"]
pipeline_LinearRegression = Pipeline([
    ('change_country_name', ChangeNameCountry()),
    ('add_regions', AddRegion()),
    ('change_distribution_2power', TransPower(['Social support', 'Freedom to make life choices'])),
    ('change_distribution_root', TransPower('Generosity', 0.5)),
    ('change_to_binary', ChangeToBinnary('Perceptions of corruption')),
    ('delete_country_col', ColumnDropper(['Country', 'Country_cc', 'Region','Regional indicator', 'Year', 'Ranking'])),
    ('one_hot_region', OneH('Region_cc')),
    ('scal_all_data', StandardScaler()),
    ('model', LinearRegression())
])

# Uruchomienie procesu treningu modelu
pipeline_LinearRegression.fit(X_train_1, y_train_1)

Pipeline(steps=[('change_country_name', ChangeNameCountry()),
                ('add_regions', AddRegion()),
                ('change_distribution_2power',
                 TransPower(column=['Social support',
                                    'Freedom to make life choices'])),
                ('change_distribution_root',
                 TransPower(column='Generosity', power=0.5)),
                ('change_to_binary',
                 ChangeToBinnary(columns='Perceptions of corruption')),
                ('delete_country_col',
                 ColumnDropper(columns=['Country', 'Country_cc', 'Region',
                                        'Regional indicator', 'Year',
                                        'Ranking'])),
                ('one_hot_region', OneH(column='Region_cc')),
                ('scal_all_data', StandardScaler()),
                ('model', LinearRegression())])

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

# Wyliczanie predykcji
y_train_pred_LinearRegression = pipeline_LinearRegression.predict(X_train_1)
y_test_pred_LinearRegression = pipeline_LinearRegression.predict(X_test_1)

# Zapisanie wyników
add_model_summary("Linear Regression", y_train_1, y_train_pred_LinearRegression, y_test_1, y_test_pred_LinearRegression)

Każdy kolejny model uczenia maszynowego wymaga dostrojenia hiperparametrów, aby osiągnąć najlepsze wyniki. W tym celu, dla wszystkich poniższych algorytmów, zastosowano metodę przeszukiwania losowego, która pozwoliła automatycznie dobrać optymalne zestawy parametrów.

Model SVM

# Tworzenie pipeline do drzewa decyzyjnego [Wykorzystano autorskie transformery zdefiniowane w zewnętrznym pliku "classes_preprocessing.py"]
pipeline_SVM= Pipeline([
    ('change_country_name', ChangeNameCountry()),
    ('add_regions', AddRegion()),
    ('change_distribution_2power', TransPower(['Social support', 'Freedom to make life choices'])),
    ('change_distribution_root', TransPower('Generosity', 0.5)),
    ('change_to_binary', ChangeToBinnary('Perceptions of corruption')),
    ('delete_country_col', ColumnDropper(['Country', 'Country_cc', 'Region', 'Regional indicator', 'Year', 'Ranking'])),
    ('one_hot_region', OneH('Region_cc')),
    ('scal_all_data', StandardScaler()),
    ('model', SVR())])

# Tworzenie siatki hiperparametrów do losowego przeszukiwania
param_grid_SVM = {
    'model__C': [0.1, 1, 10, 100],
    'model__epsilon': [0.01, 0.1, 0.5, 1],
    'model__gamma': ['scale', 'auto', 0.1, 1]}

# Konfiguracja przeszukiwania siatki
SVM = GridSearchCV(pipeline_SVM, param_grid_SVM, cv=3, scoring='r2')
# Uruchomienie procesu strojenia modelu
SVM.fit(X_train_1, y_train_1)

GridSearchCV(cv=3,
             estimator=Pipeline(steps=[('change_country_name',
                                        ChangeNameCountry()),
                                       ('add_regions', AddRegion()),
                                       ('change_distribution_2power',
                                        TransPower(column=['Social support',
                                                           'Freedom to make '
                                                           'life choices'])),
                                       ('change_distribution_root',
                                        TransPower(column='Generosity',
                                                   power=0.5)),
                                       ('change_to_binary',
                                        ChangeToBinnary(columns='Perceptions '
                                                                'of '
                                                                'corruption')),
                                       ('delete_country_col',
                                        ColumnDropper(columns=['Country',
                                                               'Country_cc',
                                                               'Region',
                                                               'Regional '
                                                               'indicator',
                                                               'Year',
                                                               'Ranking'])),
                                       ('one_hot_region',
                                        OneH(column='Region_cc')),
                                       ('scal_all_data', StandardScaler()),
                                       ('model', SVR())]),
             param_grid={'model__C': [0.1, 1, 10, 100],
                         'model__epsilon': [0.01, 0.1, 0.5, 1],
                         'model__gamma': ['scale', 'auto', 0.1, 1]},
             scoring='r2')

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

# Wyświetlenie najlepszych znalezionych hiperparametrów dla modelu
print("Najlepsze znalezione hiperparametry dla modelu SVM:")
print(SVM.best_params_)

Najlepsze znalezione hiperparametry dla modelu SVM:
{'model__C': 10, 'model__epsilon': 0.1, 'model__gamma': 'auto'}

# Wyliczanie predykcji
y_train_pred = SVM.predict(X_train_1)
y_test_pred = SVM.predict(X_test_1)

# Zapisanie wyników modelu
add_model_summary("SVM", y_train_1, y_train_pred, y_test_1, y_test_pred)

Model drzewa decyzyjnego

# Tworzenie pipeline do drzewa decyzyjnego [Wykorzystano autorskie transformery zdefiniowane w zewnętrznym pliku "classes_preprocessing.py"]
pipeline_DecisionTree= Pipeline([
    ('change_country_name', ChangeNameCountry()),
    ('add_regions', AddRegion()),
    ('change_distribution_2power', TransPower(['Social support', 'Freedom to make life choices'])),
    ('change_distribution_root', TransPower('Generosity', 0.5)),
    ('change_to_binary', ChangeToBinnary('Perceptions of corruption')),
    ('delete_country_col', ColumnDropper(['Country', 'Country_cc', 'Region', 'Regional indicator', 'Year', 'Ranking'])),
    ('one_hot_region', OneH('Region_cc')),
    ('scal_all_data', StandardScaler()),
    ('model', DecisionTreeRegressor(random_state=42))])

# Tworzenie siatki hiperparametrów do losowego przeszukiwania
param_grid_DecisionTree = {
    'model__max_depth': [None, 3, 5, 7, 10, 15, 20],
    'model__min_samples_split': [2, 5, 10, 20],
    'model__min_samples_leaf': [1, 2, 4, 6]}

# Konfiguracja przeszukiwania siatki
DecisionTree = GridSearchCV(pipeline_DecisionTree, param_grid_DecisionTree, cv=3, scoring='r2')
# Uruchomienie procesu strojenia modelu
DecisionTree.fit(X_train_1, y_train_1)

GridSearchCV(cv=3,
             estimator=Pipeline(steps=[('change_country_name',
                                        ChangeNameCountry()),
                                       ('add_regions', AddRegion()),
                                       ('change_distribution_2power',
                                        TransPower(column=['Social support',
                                                           'Freedom to make '
                                                           'life choices'])),
                                       ('change_distribution_root',
                                        TransPower(column='Generosity',
                                                   power=0.5)),
                                       ('change_to_binary',
                                        ChangeToBinnary(columns='Perceptions '
                                                                'of '
                                                                'corr...
                                        ColumnDropper(columns=['Country',
                                                               'Country_cc',
                                                               'Region',
                                                               'Regional '
                                                               'indicator',
                                                               'Year',
                                                               'Ranking'])),
                                       ('one_hot_region',
                                        OneH(column='Region_cc')),
                                       ('scal_all_data', StandardScaler()),
                                       ('model',
                                        DecisionTreeRegressor(random_state=42))]),
             param_grid={'model__max_depth': [None, 3, 5, 7, 10, 15, 20],
                         'model__min_samples_leaf': [1, 2, 4, 6],
                         'model__min_samples_split': [2, 5, 10, 20]},
             scoring='r2')

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

# Wyświetlenie najlepszych znalezionych hiperparametrów dla modelu
print("Najlepsze znalezione hiperparametry dla modelu drzewa decyzyjnego:")
print(DecisionTree.best_params_)

Najlepsze znalezione hiperparametry dla modelu drzewa decyzyjnego:
{'model__max_depth': 7, 'model__min_samples_leaf': 2, 'model__min_samples_split': 5}

# Wyliczanie predykcji
y_train_pred = DecisionTree.predict(X_train_1)
y_test_pred = DecisionTree.predict(X_test_1)
# Zapisanie wyników modelu
add_model_summary("DecisionTree", y_train_1, y_train_pred, y_test_1, y_test_pred)

Model lasu losowego

# Tworzenie pipeline do lasu losowego [Wykorzystano autorskie transformery zdefiniowane w zewnętrznym pliku "classes_preprocessing.py"]
pipeline_RandomForest = Pipeline([
    ('change_country_name', ChangeNameCountry()),
    ('add_regions', AddRegion()),
    ('change_distribution_2power', TransPower(['Social support', 'Freedom to make life choices'])),
    ('change_distribution_root', TransPower('Generosity', 0.5)),
    ('change_to_binary', ChangeToBinnary('Perceptions of corruption')),
    ('delete_country_col', ColumnDropper(['Country', 'Country_cc', 'Region', 'Regional indicator', 'Year', 'Ranking'])),
    ('One_hot_region', OneH('Region_cc')),
    ('scal_all_data', StandardScaler()),
    ('model', RandomForestRegressor(random_state=42))])

# Tworzenie siatki hiperparametrów do losowego przeszukiwania
param_grid_RandomForest = {
    'model__n_estimators': [10, 100, 150, 200, 250, 300, 400, 500, 750, 1000],
    'model__max_depth': [None, 3, 5, 7, 9, 11, 13, 14],}
# Konfiguracja przeszukiwania siatki
Random_Forest = GridSearchCV(pipeline_RandomForest, param_grid_RandomForest, cv=3, scoring='r2')

# Uruchomienie procesu strojenia modelu
Random_Forest.fit(X_train_1, y_train_1)

GridSearchCV(cv=3,
             estimator=Pipeline(steps=[('change_country_name',
                                        ChangeNameCountry()),
                                       ('add_regions', AddRegion()),
                                       ('change_distribution_2power',
                                        TransPower(column=['Social support',
                                                           'Freedom to make '
                                                           'life choices'])),
                                       ('change_distribution_root',
                                        TransPower(column='Generosity',
                                                   power=0.5)),
                                       ('change_to_binary',
                                        ChangeToBinnary(columns='Perceptions '
                                                                'of '
                                                                'corr...
                                        ColumnDropper(columns=['Country',
                                                               'Country_cc',
                                                               'Region',
                                                               'Regional '
                                                               'indicator',
                                                               'Year',
                                                               'Ranking'])),
                                       ('One_hot_region',
                                        OneH(column='Region_cc')),
                                       ('scal_all_data', StandardScaler()),
                                       ('model',
                                        RandomForestRegressor(random_state=42))]),
             param_grid={'model__max_depth': [None, 3, 5, 7, 9, 11, 13, 14],
                         'model__n_estimators': [10, 100, 150, 200, 250, 300,
                                                 400, 500, 750, 1000]},
             scoring='r2')

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

# Wyświetlenie najlepszych znalezionych hiperparametrów dla modelu
print("Najlepsze znalezione hiperparametry dla modelu lasu losowego:")
print(Random_Forest.best_params_)

Najlepsze znalezione hiperparametry dla modelu lasu losowego:
{'model__max_depth': None, 'model__n_estimators': 750}

# Wyliczanie predykcji
y_train_pred = Random_Forest.predict(X_train_1)
y_test_pred = Random_Forest.predict(X_test_1)

# Zapisanie wyników modelu
add_model_summary("Random_Forest", y_train_1, y_train_pred, y_test_1, y_test_pred)

Model XGBoost

# Tworzenie pipeline do XGBoost [Wykorzystano autorskie transformery zdefiniowane w zewnętrznym pliku "classes_preprocessing.py"]
pipeline_XGBoost = Pipeline([
    ('change_country_name', ChangeNameCountry()),
    ('add_regions', AddRegion()),
    ('change_distribution_2power', TransPower(['Social support', 'Freedom to make life choices'])),
    ('change_distribution_root', TransPower('Generosity', 0.5)),
    ('change_to_binary', ChangeToBinnary('Perceptions of corruption')),
    ('delete_country_col', ColumnDropper(['Country', 'Country_cc', 'Region', 'Regional indicator', 'Year', 'Ranking'])),
    ('One_hot_region', OneH('Region_cc')),
    ('scal_all_data', StandardScaler()),
    ('model', XGBRegressor(random_state=42))])

# Siatka hiperparametrów do losowego przeszukiwania
param_grid_XGBoost = {
    'model__n_estimators': [100, 200, 500],
    'model__max_depth': [3, 5, 7],
    'model__learning_rate': [0.01, 0.05, 0.1]}
# Konfiguracja przeszukiwania siatki
XGBoost = GridSearchCV(pipeline_XGBoost, param_grid_XGBoost, cv=3, scoring='r2')

# Uruchomienie procesu strojenia modelu
XGBoost.fit(X_train_1, y_train_1)

GridSearchCV(cv=3,
             estimator=Pipeline(steps=[('change_country_name',
                                        ChangeNameCountry()),
                                       ('add_regions', AddRegion()),
                                       ('change_distribution_2power',
                                        TransPower(column=['Social support',
                                                           'Freedom to make '
                                                           'life choices'])),
                                       ('change_distribution_root',
                                        TransPower(column='Generosity',
                                                   power=0.5)),
                                       ('change_to_binary',
                                        ChangeToBinnary(columns='Perceptions '
                                                                'of '
                                                                'corr...
                                                     max_cat_to_onehot=None,
                                                     max_delta_step=None,
                                                     max_depth=None,
                                                     max_leaves=None,
                                                     min_child_weight=None,
                                                     missing=nan,
                                                     monotone_constraints=None,
                                                     multi_strategy=None,
                                                     n_estimators=None,
                                                     n_jobs=None,
                                                     num_parallel_tree=None,
                                                     random_state=42, ...))]),
             param_grid={'model__learning_rate': [0.01, 0.05, 0.1],
                         'model__max_depth': [3, 5, 7],
                         'model__n_estimators': [100, 200, 500]},
             scoring='r2')

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

# Wyświetlenie najlepszych znalezionych hiperparametrów dla modelu
print("Najlepsze znalezione hiperparametry dla modelu XGBoost:")
print(XGBoost.best_params_)

Najlepsze znalezione hiperparametry dla modelu XGBoost:
{'model__learning_rate': 0.1, 'model__max_depth': 5, 'model__n_estimators': 500}

# Wyliczanie predykcji
y_train_pred = XGBoost.predict(X_train_1)
y_test_pred = XGBoost.predict(X_test_1)

# Zapisanie wyników modelu
add_model_summary("XGBoost", y_train_1, y_train_pred, y_test_1, y_test_pred)

Model sieci neuronowej (MLP)

# Definicja architektury modelu MLP i hiperparametrów do optymalizacji
def build_model(hp):
    n_hidden = hp.Int("n_hidden", min_value=1, max_value=6, default=1)
    n_neurons = hp.Int("n_neurons", min_value=2, max_value=256)
    learning_rate = hp.Float("learning_rate", min_value=1e-4, max_value=1e-2, sampling="log")
    optimizer = hp.Choice("optimizer", values=["sgd", "adam"])
    if optimizer == "sgd":
        optimizer = tf.keras.optimizers.SGD(learning_rate=learning_rate)
    else:
        optimizer = tf.keras.optimizers.Adam(learning_rate=learning_rate)
    
    model = tf.keras.Sequential()
    model.add(tf.keras.layers.Flatten())
    for _ in range(n_hidden):
        model.add(tf.keras.layers.Dense(n_neurons, activation="relu"))
    model.add(tf.keras.layers.Dense(1))
    model.compile(loss="mse", optimizer=optimizer, metrics=["mae"])
    return model

# Tworzenie pipeline dla sieci nerunowej [Wykorzystano autorskie transformery zdefiniowane w zewnętrznym pliku "classes_preprocessing.py"]
preprocessing_MLP = Pipeline([
    ('change_country_name', ChangeNameCountry()),
    ('add_regions', AddRegion()),
    ('change_distribution_2power', TransPower(['Social support', 'Freedom to make life choices'])),
    ('change_distribution_root', TransPower('Generosity', 0.5)),
    ('change_to_binary', ChangeToBinnary('Perceptions of corruption')),
    ('delete_country_col', ColumnDropper(['Country', 'Country_cc', 'Region', 'Regional indicator', 'Year', 'Ranking'])),
    ('One_hot_region', OneH('Region_cc')),
    ('scal_all_data', StandardScaler()),])

# Preprocessing danych treningowych i testowych
X_train_1_pre = preprocessing_MLP.fit_transform(X_train_1)
X_test_1_pre = preprocessing_MLP.transform(X_test_1)

# Definicja funkcji tworzącej nazwy folderów z datą, aby nie nadpisywać starych wyników
def get_run_logdir(root_logdir="logs_MLP"):
    return Path(root_logdir) / strftime("run_%Y_%m_%d_%H_%M_%S")
# Konfiguracja zapisu logów do wizualizacji w TensorBoard
run_logdir = get_run_logdir()
tensorboard_cb = tf.keras.callbacks.TensorBoard(run_logdir)

# Konfiguracja przeszukiwania
random_search_tuner = kt.RandomSearch(
    build_model, objective="val_mae", max_trials=5, overwrite=True,
    directory="project_happines", project_name="happines_neuron_random_search", seed=42)
# Uruchomienie procesu strojenia modelu
random_search_tuner.search(X_train_1_pre, y_train_1, epochs=10, validation_split=0.2, callbacks=[tensorboard_cb])

Trial 5 Complete [00h 00m 07s]
val_mae: 0.8084776997566223

Best val_mae So Far: 0.31277671456336975
Total elapsed time: 00h 00m 47s

# Wyświetlenie najlepszych znalezionych hiperparametrów dla modelu MLP
print("Najlepsze znalezione hiperparametry dla modelu MLP:")
print(random_search_tuner.get_best_hyperparameters()[0].values)

Najlepsze znalezione hiperparametry dla modelu MLP:
{'n_hidden': 3, 'n_neurons': 64, 'learning_rate': 0.00905127409782462, 'optimizer': 'adam'}

Na poniższym wykresie przedstawiono historię uczenia najlepszego wariantu sieci neuronowej, wyłonionego w procesie losowego przeszukiwania. Linia pomarańczowa przedstawia funkcję straty (MSE) dla danych treningowych, natomiast linia niebieska odpowiada danym walidacyjnym. Proces uczenia przebiegł prawidłowo, po gwałtownym spadku błędu w pierwszej epoce dla danych treningowych, już od 3 epoki obie krzywe zbiegają się i stabilizują na bardzo niskim poziomie. Brak istotnych wahań oraz fakt, że w końcowej fazie obie linie niemal całkowicie się pokrywają, świadczy o braku zjawisku przeuczenia czy niedotrenowania.

# Wyciszenie ostrzeżeń dla tej komórki (komunikaty nie wpływają na wyniki, a obniżają czytelność)
with warnings.catch_warnings():
    warnings.simplefilter("ignore")
    # Wybieranie najlepszgo modelu
    best_mlp_model = random_search_tuner.get_best_models()[0]
    # Wyliczanie predykcji
    y_train_pred_MLP = best_mlp_model.predict(X_train_1_pre)
    y_test_pred_MLP = best_mlp_model.predict(X_test_1_pre)
    # Zapisanie wyników modelu
    add_model_summary("MLP", y_train_1, y_train_pred_MLP, y_test_1, y_test_pred_MLP)

22/22 ━━━━━━━━━━━━━━━━━━━━ 0s 8ms/step
6/6 ━━━━━━━━━━━━━━━━━━━━ 0s 4ms/step 

7. Porównanie i omówienie wyników modeli

# Tworzenie ramki danych z listy wyników
df_results = pd.DataFrame(model_results)
# Ustawienie kolumny model jako indeks ramki
df_results.set_index("Model", inplace=True)
# Wyświetlenie ramki danych
display(df_results.sort_values(by="Test R2", ascending=False).round(4))

	Train R2	Test R2	Train RMSE	Test RMSE	Train MAE	Test MAE
Model
XGBoost	0.9998	0.9177	0.0166	0.3631	0.0116	0.2684
Random_Forest	0.9854	0.9172	0.1299	0.3642	0.0976	0.2747
SVM	0.9346	0.9012	0.2751	0.3978	0.1854	0.2852
Linear Regression	0.8490	0.8619	0.4181	0.4703	0.3195	0.3504
MLP	0.8833	0.8493	0.3676	0.4914	0.2773	0.3689
DecisionTree	0.9063	0.8384	0.3294	0.5088	0.2347	0.3825

Na podstawie powyższej tabeli można wyciągnąc następujące wnioski:

• Najlepszy model - najlepszym modelem okazał się XGBoost, który na zbiorze testowym wyjaśnia ponad 91% wariancji. Średni błąd bezwzględny (MAE) dla danych testowych w przybliżeniu wyniósł zaledwie 0.27 punktu. Oznacza to, że na 10-stopniowej skali szczęścia nasze predykcje mylą się średnio o mniej niż 0.3, co jest wynikiem bardzo zadowalającym. Wyniki te pozwalają twierdząco odpowiedzieć na postawione we wstępie pytanie badawcze. Okazuje się, że wykorzystując zaledwie 7 zmiennych, jesteśmy w stanie z wysoką dokładnością przewidzieć, jak mieszkańcy danego kraju ocenią swój poziom szczęścia.
• Stabilność i Efektywność Regresji Liniowej: - model Regresji Liniowej osiągnął bardzo wysoki i przede wszystkim wyjątkowo stabilny wynik wyjaśniając około 85% wariancji zarówno dla zbioru testowego i treningowego. Tak wysoka skuteczność najprostszego algorytmu dowodzi, że zależności między zmiennymi a poziomem szczęścia ma w dużej mierze strukturę liniową (co jest zgodne z wnioskami ze zbudowanej wcześniej macierzy korelacji Pearsona).
• Złożoność a efektywność - Słabszy wynik od regresji liniowej w zestawieniu osiągnęła sieć neuronowa (MLP). Zbiór danych charakteryzował się niską wymiarowością (zaledwie 6 zmiennych, w tym 5 numerycznych) oraz silną strukturą liniowych zależności. Sieci neuronowe zazwyczaj pokazują swoją przewagę w złożonych, nieliniowych problemach przy dużej ilości danych. Wyniki potwierdzają regułę, że w przypadku prostszych, liniowych problemów, klasyczne modele uczenia maszynowego są rozwiązaniem bardziej efektywnym i stabilnym, przy czym są także prostsze do wdrożenia.

8. Analiza SHAP dla najlepszego modelu

# Wybieranie pipeline najlepszego modelu
pipeline_best_params = XGBoost.best_estimator_
# Wyodrębnienie transformatora danych
preprocessor = pipeline_best_params[:-1] 
# Wyodrębnienie modelu
model = pipeline_best_params[-1] 
# Preprocessing danych
X_train_preprocessed = preprocessor.transform(X_train_1)
# Pobranie nazw kolumn
cols = preprocessor[-1].get_feature_names_out()
# Tworzenie explainera danych
explainer = shap.TreeExplainer(model)
# Obliczenie wartości SHAP
shap_values = explainer.shap_values(X_train_preprocessed)
# Wyświetlenie wykresu
shap.summary_plot(shap_values, X_train_preprocessed, feature_names=cols)

W ramach analizy stworzono także wykres wartości SHAP. Wizualizacje nalezy intepretować w następujący sposób:

• Oś Y (Ranking istotności): Cechy są uszeregowane pionowo według ich całkowitego wpływu na model. Zmienne znajdujące się na samej górze mają statystycznie największe znaczenie dla przewidywanego poziomu szczęścia.

• Oś X (Kierunek i siła wpływu): Każdy punkt na wykresie reprezentuje pojedynczy kraj. Punkty po lewej stronie osi zero (ujemne wartości SHAP) oznaczają, że dana cecha obniża prognozowany poziom szczęścia, natomiast punkty po prawej (dodatnie wartości SHAP) – że go podwyższają. Kluczowa jest tutaj odległość od środka: im dalej w lewo znajduje się punkt, tym silniejszy jest jego negatywny wpływ na wynik, i analogicznie – im wyższa dodatnia wartość SHAP (im dalej w prawo), tym mocniej dana cecha podbija prognozę szczęścia.

• Kolor (Wartość zmiennej): Skala kolorystyczna pozwala zrozumieć relację między wartością cechy a jej efektem. Kolor czerwony reprezentuje wysokie wartości danej zmiennej, a niebieski wartości niskie.

Na podstawie powyższego wykresu można zweryfikować postawione na wstępie hipotezy:

• Zamożność społeczeństwa, wyrażona przez PKB na mieszkańca skorygowana o parytet siły nabywczej, stanowi jeden z najsilniejszy predyktorów poziomu szczęścia:

  Jak widać na podstawie powyższego wykresu SHAP zmienna GDP per capita (Produkt krajowy brutto na mieszkańca skorygowany o parytet siły nabywczej) ma największy wpływ na wartość predykcji, co potwierdza hipotezę jakoby zmienna ta z dostępnego zbioru zmiennych była jednym z najlepszym predyktorem szczęścia.

• Hojność jest najsłabszym predyktorem szczęścia:

  Analiza wykazała, że hipoteza ta nie znajduje potwierdzenia w danych. Zmienna Generosity (Hojność) uplasowała się w środkowej części rankingu, wykazując silniejszy wpływ na wyniki modelu niż między innymi postrzeganie korupcji (Perceptions of corruption) czy oczekiwana długość życia w zdrowiu (Healthy life expectancy). Oznacza to, że hojność społeczeństwa ma realne i zauważalne znaczenie dla modelu predykcyjnego i zdecydowanie nie jest najgorszym predyktorem z dostępnego zbioru

Ciekawymi obserwacjami są także:

• Istotna rola przynależności regionalnej w predykcji: Analiza SHAP potwierdziła wnioski z wcześniejszej eksploracyjnej analizy danych. Zaobserwowano, że regiony sklasyfikowane na wykresach pudełkowych (boxplot) jako obszary o skrajnym poziomie dobrostanu (bardzo wysokim lub bardzo niskim, z wyłączeniem środkowej i zachodniej Afryki), mają decydujący wpływ na wynik predykcji. W przypadku krajów należących do regionów o przeciętnym poziomie szczęścia (oraz Afryki środkowej i zachodniej), wpływ zmiennej regionalnej na model jest znikomy.
• Korupcja jako najsłabszy predyktor: Zmienna Perceptions of corruption (postrzeganie korupcji) okazała się najmniej wpływową cechą w modelu. Skupienie wszystkich punktów danych wokół wartości zero na osi SHAP dowodzi, że czynnik ten ma marginalne znaczenie dla końcowej oceny szczęścia w porównaniu do pozostałych zmiennych.
• Zaskakująco niewielki wpływ śrdniej długości życia w zdrowiu: Zmienna Healthy life expectancy (oczekiwana długośc życia w zdrowiu) wykazuje słabszy wpływ na predykcję niż Freedom to make life choices (Postrzegana swoboda podejmowania ważnych decyzji życiowych) czy Generosity (hojność). Jest to nieoczywista obserwacja i pokazuje istotność czynników społecznych w ocenie szczęścia.

9. Aplikacja

#Zapisanie najlepszego modelu aby wykorzystać go w aplikacji
joblib.dump(XGBoost, 'XGBoost_happines')

['XGBoost_happines']

W ramach projektu stworzono także aplikację będącą praktyczną implementacją modelu, która służy do interaktywnej weryfikacji postawionych hipotez oraz testowania zachowania algorytmu dla dowolnych danych wejściowych.
Kliknij, aby przestestować aplikację

10. Podsumowanie

Osiągnięto główny cel projektu: Zbudowany model uczenia maszynowego osiągnął bardzo satysfakcjonujący wynik, wyjaśniając ponad 91% wariancji na zbiorze testowym. Świadczy to o wysokiej dokładności modelu w predykcjach.

Weryfikacja hipotez i pytań badawczych: W toku prac udało się odpowiedzieć na wszystkie postawione pytania oraz zweryfikować hipotezy w następujący sposób:

• Zamożność społeczeństwa, wyrażona przez PKB na mieszkańca skorygowana o parytet siły nabywczej, stanowi jeden z najsilniejszy predyktorów poziomu szczęścia:

  Za pomocą analizy wartości SHAP potwierdzono, że zmienna GDP per capita (odpowiadająca średniej wartości PKB na osobę skorygowanej o parytet siły nabywczej) jest najważniejszym czynnikiem wpływającym na poziom szczęścia. Jest to więc jeden z najsilniejszych predyktorów w modelu.

• Hojność jest najsłabszym predyktorem szczęścia:

  Za pomocą analizy SHAP udało się obalić hipotezę jakoby zmienna Generosity (hojność) była najsłabszym predyktorem modelu. Pokazano, że chociażby zmienne takie jak Healthy life expectation (oczekiwana długość życia), czy Perceptions of corruption (postrzeganie korupcji) mają zdecydowany mniejszy wpływ na wynik predykcji.

• Czy pandemia COVID-19 negatywnie wpłynęła na poziom sczęścia?:

  Podczas eksploracyjnej analizy danych wykazano, że na poziomie globalnym nie wystąpiły znaczne różnice w średniej wartości szczęścia na przestrzeni badanych lat. Analiza trendu czasowego ujawniła, że nawet w okresie największych obostrzeń związanych z pandemią COVID-19 (lata 2020–2022), wskaźnik szczęścia odnotował nieznaczny wzrost. A statystyczny test ANOVA potwierdził brak istotnyh różnic między latami. Sugeruje to, że pandemia nie miała destrukcyjnego wpływu na globalne postrzeganie szczęścia u ludzi.

• Czy da się na podstawie zaledwie kilku zmiennych określić poziom szczęścia w danym państwie?:

  Tak. Osiągnięty wynik współczynnika determinacji na zbiorze testowym (0,91) jednoznacznie potwierdza, że niewielka liczba kluczowych wskaźników jest wystarczająca do zbudowania wysoce skutecznego modelu predykcyjnego. Udowodniono, że nie potrzeba wielu złożonych zmiennych. Wystarczy kilka kluczowych informacji, by bardzo precyzyjnie przewidzieć, jak szczęśliwi są mieszkańcy danego państwa.

Autor:

Jan Walkiewicz